home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / STATRAND.ZIP / BDTR.FOR next >
Text File  |  1985-12-31  |  6KB  |  246 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE BDTR
  5. C
  6. C     PURPOSE
  7. C        COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U,
  8. C        DISTRIBUTED ACCORDING TO THE BETA DISTRIBUTION WITH
  9. C        PARAMETERS A AND B, IS LESS THAN OR EQUAL TO X.  F(A,B,X),
  10. C        THE ORDINATE OF THE BETA DENSITY AT X, IS ALSO COMPUTED.
  11. C
  12. C     USAGE
  13. C        CALL BDTR(X,A,B,P,D,IER)
  14. C
  15. C     DESCRIPTION OF PARAMETERS
  16. C        X    - INPUT SCALAR FOR WHICH P(X) IS COMPUTED.
  17. C        A    - BETA DISTRIBUTION PARAMETER (CONTINUOUS).
  18. C        B    - BETA DISTRIBUTION PARAMETER (CONTINUOUS).
  19. C        P    - OUTPUT PROBABILITY.
  20. C        D    - OUTPUT DENSITY.
  21. C        IER - RESULTANT ERROR CODE WHERE
  22. C        IER= 0 --- NO ERROR
  23. C        IER=-1,+1  CDTR HAS BEEN CALLED AND AN ERROR HAS
  24. C               OCCURRED.  SEE CDTR.
  25. C        IER=-2 --- AN INPUT PARAMETER IS INVALID.  X IS LESS
  26. C               THAN 0.0 OR GREATER THAN 1.0, OR EITHER A OR
  27. C               B IS LESS THAN 0.5 OR GREATER THAN 10**(+5).
  28. C               P AND D ARE SET TO -1.E75.
  29. C        IER=+2 --- INVALID OUTPUT.  P IS LESS THAN ZERO OR
  30. C               GREATER THAN ONE.  P IS SET TO 1.E75.
  31. C
  32. C     REMARKS
  33. C        SEE MATHEMATICAL DESCRIPTION.
  34. C
  35. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  36. C        DLGAM
  37. C        NDTR
  38. C        CDTR
  39. C
  40. C     METHOD
  41. C        REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL
  42. C        DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE,
  43. C        IBM RESEARCH REPORT RC-1094, 1963.
  44. C
  45. C     ..................................................................
  46. C
  47.       SUBROUTINE BDTR(X,A,B,P,D,IER)
  48.       DOUBLE PRECISION XX,DLXX,DL1X,AA,BB,G1,G2,G3,G4,DD,PP,XO,FF,FN,
  49.      1XI,SS,CC,RR,DLBETA
  50. C
  51. C     TEST FOR VALID INPUT DATA
  52. C
  53.       IF(A-(.5-1.E-5)) 640,10,10
  54.    10 IF(B-(.5-1.E-5)) 640,20,20
  55.    20 IF(A-1.E+5) 30,30,640
  56.    30 IF(B-1.E+5) 40,40,640
  57.    40 IF(X) 640,50,50
  58.    50 IF(1.-X) 640,60,60
  59. C
  60. C     COMPUTE LOG(BETA(A,B))
  61. C
  62.    60 AA=DBLE(A)
  63.       BB=DBLE(B)
  64.       CALL DLGAM(AA,G1,IOK)
  65.       CALL DLGAM(BB,G2,IOK)
  66.       CALL DLGAM(AA+BB,G3,IOK)
  67.       DLBETA=G1+G2-G3
  68. C
  69. C     TEST FOR X NEAR 0.0 OR 1.0
  70. C
  71.       IF(X-1.E-8) 80,80,70
  72.    70 IF((1.-X)-1.E-8) 130,130,140
  73.    80 P=0.0
  74.       IF(A-1.) 90,100,120
  75.    90 D=1.E+38
  76.       GO TO 660
  77.   100 DD=-DLBETA
  78.       IF(DD+1.68D02)  120,120,110
  79.   110 DD=DEXP(DD)
  80.       D=SNGL(DD)
  81.       GO TO 660
  82.   120 D=0.0
  83.       GO TO 660
  84.   130 P=1.0
  85.       IF(B-1.) 90,100,120
  86. C
  87. C     SET PROGRAM PARAMETERS
  88. C
  89.   140 XX=DBLE(X)
  90.       DLXX=DLOG(XX)
  91.       DL1X=DLOG(1.D0-XX)
  92.       XO=XX/(1.D0-XX)
  93.       ID=0
  94. C
  95. C     COMPUTE ORDINATE
  96. C
  97.       DD=(AA-1.D0)*DLXX+(BB-1.D0)*DL1X-DLBETA
  98.       IF(DD-1.68D02) 150,150,160
  99.   150 IF(DD+1.68D02) 170,170,180
  100.   160 D=1.E38
  101.       GO TO 190
  102.   170 D=0.0
  103.       GO TO 190
  104.   180 DD=DEXP(DD)
  105.       D=SNGL(DD)
  106. C
  107. C     A OR B OR BOTH WITHIN 1.E-8 OF 1.0
  108. C
  109.   190 IF(ABS(A-1.)-1.E-8)  200,200,210
  110.   200 IF(ABS(B-1.)-1.E-8)  220,220,230
  111.   210 IF(ABS(B-1.)-1.E-8)  260,260,290
  112.   220 P=X
  113.       GO TO 660
  114.   230 PP=BB*DL1X
  115.       IF(PP+1.68D02) 240,240,250
  116.   240 P=1.0
  117.       GO TO 660
  118.   250 PP=DEXP(PP)
  119.       PP=1.D0-PP
  120.       P=SNGL(PP)
  121.       GO TO 600
  122.   260 PP=AA*DLXX
  123.       IF(PP+1.68D02) 270,270,280
  124.   270 P=0.0
  125.       GO TO 660
  126.   280 PP=DEXP(PP)
  127.       P=SNGL(PP)
  128.       GO TO 600
  129. C
  130. C     TEST FOR A OR B GREATER THAN 1000.0
  131. C
  132.   290 IF(A-1000.) 300,300,310
  133.   300 IF(B-1000.) 330,330,320
  134.   310 XX=2.D0*AA/XO
  135.       XS=SNGL(XX)
  136.       AA=2.D0*BB
  137.       DF=SNGL(AA)
  138.       CALL CDTR(XS,DF,P,DUMMY,IER)
  139.       P=1.0-P
  140.       GO TO 670
  141.   320 XX=2.D0*BB*XO
  142.       XS=SNGL(XX)
  143.       AA=2.D0*AA
  144.       DF=SNGL(AA)
  145.       CALL CDTR(XS,DF,P,DUMMY,IER)
  146.       GO TO 670
  147. C
  148. C     SELECT PARAMETERS FOR CONTINUED FRACTION COMPUTATION
  149. C
  150.   330 IF(X-.5) 340,340,380
  151.   340 IF(AA-1.D0) 350,350,360
  152.   350 RR=AA+1.D0
  153.       GO TO 370
  154.   360 RR=AA
  155.   370 DD=DLXX/5.D0
  156.       DD=DEXP(DD)
  157.       DD=(RR-1.D0)-(RR+BB-1.D0)*XX*DD +2.D0
  158.       IF(DD) 420,420,430
  159.   380 IF(BB-1.D0) 390,390,400
  160.   390 RR=BB+1.D0
  161.       GO TO 410
  162.   400 RR=BB
  163.   410 DD=DL1X/5.D0
  164.       DD=DEXP(DD)
  165.       DD=(RR-1.D0)-(AA+RR-1.D0)*(1.D0-XX)*DD +2.D0
  166.       IF(DD) 430,430,420
  167.   420 ID=1
  168.       FF=DL1X
  169.       DL1X=DLXX
  170.       DLXX=FF
  171.       XO=1.D0/XO
  172.       FF=AA
  173.       AA=BB
  174.       BB=FF
  175.       G2=G1
  176. C
  177. C     TEST FOR A LESS THAN 1.0
  178. C
  179.   430 FF=0.D0
  180.       IF(AA-1.D0) 440,440,470
  181.   440 CALL DLGAM(AA+1.D0,G4,IOK)
  182.       DD=AA*DLXX+BB*DL1X+G3-G2-G4
  183.       IF(DD+1.68D02) 460,460,450
  184.   450 FF=FF+DEXP(DD)
  185.   460 AA=AA+1.D0
  186. C
  187. C     COMPUTE P USING CONTINUED FRACTION EXPANSION
  188. C
  189.   470 FN=AA+BB-1.D0
  190.       RR=AA-1.D0
  191.       II=80
  192.       XI=DFLOAT(II)
  193.       SS=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI))
  194.       SS=SS*XO
  195.       DO 480 I=1,79
  196.       II=80-I
  197.       XI=DFLOAT(II)
  198.       DD=(XI*(FN+XI))/((RR+2.D0*XI+1.D0)*(RR+2.D0*XI))
  199.       DD=DD*XO
  200.       CC=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI))
  201.       CC=CC*XO
  202.       SS=CC/(1.D0+DD/(1.D0-SS))
  203.   480 CONTINUE
  204.       SS=1.D0/(1.D0-SS)
  205.       IF(SS) 650,650,490
  206.   490 CALL DLGAM(AA+BB,G1,IOK)
  207.       CALL DLGAM(AA+1.D0,G4,IOK)
  208.       CC=G1-G2-G4+AA*DLXX+(BB-1.D0)*DL1X
  209.       PP=CC+DLOG(SS)
  210.       IF(PP+1.68D02) 500,500,510
  211.   500 PP=FF
  212.       GO TO 520
  213.   510 PP=DEXP(PP)+FF
  214.   520 IF(ID) 540,540,530
  215.   530 PP=1.D0-PP
  216.   540 P=SNGL(PP)
  217. C
  218. C     SET ERROR INDICATOR
  219. C
  220.       IF(P) 550,570,570
  221.   550 IF(ABS(P)-1.E-7) 560,560,650
  222.   560 P=0.0
  223.       GO TO 660
  224.   570 IF(1.-P) 580,600,600
  225.   580 IF(ABS(1.-P)-1.E-7) 590,590,650
  226.   590 P=1.0
  227.       GO TO 660
  228.   600 IF(P-1.E-8) 610,610,620
  229.   610 P=0.0
  230.       GO TO 660
  231.   620 IF((1.0-P)-1.E-8) 630,630,660
  232.   630 P=1.0
  233.       GO TO 660
  234.   640 IER=-2
  235.       D=-1.E38
  236.       P=-1.E38
  237.       GO TO 670
  238.   650 IER=+2
  239.       P= 1.E38
  240.       GO TO 670
  241.   660 IER=0
  242.   670 RETURN
  243.       END
  244. 60
  245.   570 IF(1.-P) 580,600,600
  246.